// cm - translate c(M) for any overdensity value.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "zbrent.h"
#include "cm.h"

#define N 1024

double min[N], cin[N];
double mout[N], cout[N];

double rdelta,c;

extern double delta,deltal;

const double rho_crit=2.7755e11;

int readcm(void) {

	char line[80];
	int count;
	
	count=0;
	while (fgets(line,80,stdin))
	{	
		if (line[0] != '#' && strlen(line) > 1) 
		{
			sscanf(line, "%lg %lg", &min[count], &cin[count]); 
			count++;
		}
	}

	return count;
}

void writecm(double delta, int _n) 
{
	int i;
	fprintf(stdout, "# m%lg c%lg\n",delta, delta);
	for (i=0;i<_n;i++)
	{
		fprintf(stdout,"%lg %lg\n", mout[i], cout[i]);
	}
}

static double deltac(double _c, double delta)
{
    return (delta/3.0)*(pow(_c,3.0))/(log(1+_c)-_c/(1+_c));
}

// from rho(delta,c) compute r_deltal findding the root of <rho>=deltal*rho_crit

static double rdeltal(double x, double c)
{

	return -3.0*deltac(c,delta)*pow(rdelta,3.0)*(x*c+log(rdelta)*rdelta+log(rdelta)*x*c-log(rdelta+x*c)*rdelta-log(rdelta+x*c)*x*c)/pow(c,3.0)/(rdelta+x*c)/pow(x,3.0)-deltal;
}


int convertcm(void)
{

	int i,n;
	double _rdeltal,rs;
 
	fprintf(stderr,"Input overdensity = %g\n",delta);
	fprintf(stderr,"Output overdensity = %g\n",deltal);

	n=readcm();

	for (i=0;i<n;i++)
	{
		c=cin[i];
		rdelta=pow(3.0*min[i]/(4.0*M_PI*delta*rho_crit),1.0/3.0);

		rs=rdelta/c;

		_rdeltal=zbrent(rdeltal,0.01,5.0,0.00001);

		cout[i]=_rdeltal/rs;
		mout[i]=4.0/3.0*M_PI*pow(_rdeltal,3.0)*deltal*rho_crit;
	}
	
	writecm(deltal,n);

	return 0;
}
